from math import *
import matplotlib.pyplot as plt
from matplotlib import font_manager as fm
from matplotlib import rcParams
from matplotlib.ticker import MultipleLocator


def XY_Axis(x_start, x_end, y_start, y_end):
    plt.rcParams['figure.dpi'] = 1200

    fig = plt.figure(figsize=(5, 4))
    ax = fig.add_subplot(111)

    xmajor = MultipleLocator(0.5 * x_end)
    ymajor_1 = MultipleLocator(y_end * 0.25)

    ax.xaxis.set_major_locator(xmajor)
    ax.tick_params(axis="x", direction="out", which="major", labelsize=9, length=3)
    plt.xticks([x_start, (x_start + x_end)/2, x_end], ['0', 'arcsin($r/L$)', '2arcsin($r/L$)'])

    ax.set_yticks([])

    plt.text(x_end * 0.488, y_start - 0.57, "θ", fontproperties = fm.FontProperties(fname = 'C:/Windows/Fonts/TIMESI.TTF', size = 14, style = 'italic'))
    plt.text(x_start - 0.09, y_end * 0.465, "f$_{3}$(θ)", fontproperties = fm.FontProperties(fname = 'C:/Windows/Fonts/TIMESI.TTF', size = 14, style = 'italic'), rotation = 90)


def D_PDF_OneWP(L, r, mAngle, N, Ni):
    xArray = []
    yArray = []
    L2 = L * L
    r2 = r * r
    for i in range(-N, 1):
        a = mAngle * (i / N)
        t = tan(a)
        b = L * tan(a)
        b2 = b * b
        s = sin(a)
        c = cos(a)
        bc = b * c
        c2 = c * c
        Esqrt = sqrt(r2 - b2 * c2)
        coorX = c * Esqrt - s * bc
        if abs(coorX) > 0.0:
            coorY = s * Esqrt + b * c2
            Fract = (coorY / coorX)
            if coorX > 0.0:
                B0 = atan(Fract)
            else:
                B0 = pi + atan(Fract)
        else:
            B0 = pi * 1.5
                
        DeltaB = acos(abs(bc / r)) * 2
        InterL = DeltaB / Ni
        B1 = B0 - DeltaB
        SumofFunct = 0.0
        for j in range(0, Ni):
            B = B1 + DeltaB * ((j + 0.5) / Ni)
            Tb = b / (sin(B) - t * cos(B))
            SumofFunct += pow(r2 - Tb * Tb, 3)
        y = (SumofFunct * InterL) / ( 2 * pi * pow(r2, 3))
        xArray.append(a + mAngle)
        yArray.append(y)
    for i in range(1, N + 1):
        a = mAngle * (i / N)
        t = tan(a)
        b = L * tan(a)
        b2 = b * b
        s = sin(a)
        c = cos(a)
        bc = b * c
        c2 = c * c
        Esqrt = sqrt(r2 - b2 * c2)
        coorX = c * Esqrt - s * bc
        if abs(coorX) > 0.0:
            coorY = s * Esqrt + b * c2
            Fract = (coorY / coorX)
            if coorX > 0.0:
                B0 = atan(Fract)
            else:
                B0 = pi - atan(-Fract)
        else:
            B0 = pi * 0.5
                
        DeltaB = acos(abs(bc / r)) * 2
        InterL = DeltaB / Ni
        B1 = B0
        SumofFunct = 0.0
        for j in range(0, Ni):
            B = B1 + DeltaB * ((j + 0.5) / Ni)
            Tb = b / (sin(B) - t * cos(B))
            SumofFunct += pow(r2 - Tb * Tb, 3)
        y = (SumofFunct * InterL) / ( 2 * pi * pow(r2, 3))
        xArray.append(a + mAngle)
        yArray.append(1.0 - y)

    xx, yy = [], []
    InterlX = mAngle / N
    LxyArray = len(xArray)
    YYmax = 0.0
    for i in range(1, LxyArray - 1):
        xx.append(xArray[i])
        Tyy = 0.5 * (yArray[i+1] - yArray[i-1]) / InterlX
        yy.append(Tyy)
        if YYmax < Tyy:
            YYmax = Tyy
    XY_Axis(0.0, 2 * mAngle, 0.0, YYmax)
    plt.plot(xx, yy, "k-", linewidth=1, linestyle="-")
    plt.savefig("Fig5.jpg")
    plt.show()


L = 100
r = 35
mAngle = asin(r / L)
N = 200
Ni = 1000
D_PDF_OneWP(L, r, mAngle, N, Ni)
